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Abstract 

A  simulation  of  a  triaxial  shaker  system  is  being  con¬ 
structed  to  evaluate  system  capabilities  and  limitations, 
supplement  live  vibration  testing,  support  test  method¬ 
ology  and  test  procedure  development,  and  aid  in  making 
controller  software  modifications.  Two  different  types  of 
models  have  been  implemented  in  the  simulation.  The 
first  is  a  physically-based  model  derived  from  a  finite  el¬ 
ement  analysis  together  with  a  model-updating  system 
identification  scheme;  the  second  is  a  parametric  model 
without  direct  physical  significance.  The  advantages  and 
disadvantages  of  each  model  for  this  application  are  con¬ 
sidered.  The  methods  used  to  derive  each  of  the  models 
are  described.  Test  cases  were  used  to  evaluate  the  effec¬ 
tiveness  of  some  of  the  methods.  Results  of  the  system 
identification  process  are  discussed.  Certain  methods  are 
found  to  produce  models  that  are  in  good  agreement  with 
measured  response  data  from  the  actual  shaker  system. 

1.  Introduction 

The  United  States  Air  Force  has  developed  a  triaxial 
shaker  system  for  shock  and  vibration  testing^.  The 
system  consists  of  a  table  upon  which  the  test  article 
is  mounted  and  eight  actuators  placed  along  three  or¬ 
thogonal  axes  (see  figure  1).  The  shaker  system  is  de¬ 
signed  to  produce  vibrations  ranging  in  frequency  from 
1  to  2500  Hz.  The  controller  for  the  system  is  capa¬ 
ble  of  multiple  modes  of  operation,  including  sine  wave, 
transient  waveform,  random  waveform,  random  on  ran¬ 
dom,  sine  on  random,  and  mission  siniulation.  A  mul¬ 
tiaxis  vibration  test  facility  has  the  potential  to  provide 
great  flexibility  in  test  design.  In  particular,  the  capa¬ 
bility  to  impart  both  translational  and  rotational  motion 
on  an  article  is  a  significant  advance  in  test  technology. 
At  the  same  time,  active  control  of  this  multiple- input, 
multiple-output  system  represents  a  formidable  problem. 

An  engineering  simulation  of  the  triaxial  shaker  sys¬ 
tem  is  being  constructed  in  order  to  evaluate  system 
capabilities  and  limitations,  supplement  live  vibration 
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testing,  support  test  methodology  and  test  procedure 
development,  and  aid  in  making  controller  software 
modifications  .  The  simulation  is  described  in  section  2. 
Two  different  types  of  models  have  been  implemented 
in  the  simulation.  The  first  is  a  physically-based  model 
derived  from  a  finite  element  analysis  together  with  a 
model-updating  system  identification  scheme^  the 
second  is  a  parametric  model  without  direct  physical 
significance^ .  The  advantages  and  disadvantages  of  each 
model  for  this  application  are  considered  here. 

System  identification  schemes  were  employed  to  ensure 
that  the  models  used  in  the  simulation  were  predictive 
of  actual  test  data.  For  the  physically-based  models, 
model-updating  system  identification  schemes  were  used. 
Various  approaches  to  model  updating  have  been  pro¬ 
posed.  For  instance,  Minas  and  Inman  systematically 
adjusted  a  state  space  model  based  on  measured  resonant 
frequencies  and  damping  ratios  (modal  data)  by  using  a 
pole  placement  method^.  Numerous  other  techniques 
have  been  cited  in  the  literature.  The  chief  advantage 
of  this  class  of  approaches  is  that  a  priori  knowledge  of 
the  system  can  be  incorporated  into  the  identification 
process. 

Many  approaches  to  updating  models  can  be  classified 
into  one  of  two  groups,  input  (or  equation)  error  for¬ 
mulations  and  output  error  formulations^  ^  Input 
error  formulations  produce  a  system  of  equations  that 
may  be  linear  in  the  parameters  to  be  identified,  if  mea¬ 
sured  responses  corresponding,  to  each  degree  of  freedom 
in  the  model  are  available  and  if  the  system  equations 
are  linear  in  the  parameters.  Output  error  formulations 
produce  a  system  of  equations  that  are  nonlinear  in  the 
parameters  to  be  identified,  but  they  can  easily  handle 
cases  were  there  are  fewer  measurements  than  degrees  of 
freedom  in  the  model  (the  usual  case).  A  difficulty  with 
the  output  error  formulation  is  successfully  converging 
on  the  true  values  of  the  parameters,  since  the  recep- 
tance  function  (which  appears  in  the  formulation)  can 
be  a  highly  nonlinear  or  even  discontinuous  function  of 
the  parameters^ 

Much  work  has  been  done  using  input  error  formula- 
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tions.  Most  updating  methods  associated  with  these 
formulations  require  one-to-one  correspondence  between 
the  analytical  degrees  of  freedom  and  the  experimental 

degrees  of  freedom  (the  measurements).  This  was  the 

11  19 

case  considered  by  Mottershead  ,  Mottershead  et  al  , 

and  Priswell  .  In  cases  where  there  are  fewer  measure¬ 
ments  than  analytical  degrees  of  freedom,  the  one-to-one 
correspondence  can  be  accomplished  through  reduction 
of  the  analytical  system  matrices  or  expansion  of  the  ex¬ 
perimental  data.  Methods  for  dynamic  reduction  were 
presented  by  Lammens  et  al^^  and  by  D’Amabrogio  and 
Fregolent^^.  These  methods  assume  that  no  forces  are 
applied  at  unmeasured  degrees  of  freedom.  This  assump¬ 
tion  is  not  valid  for  the  triaxial  shaker  problem;  an  ap¬ 
proach  that  is  valid  for  the  case  of  the  triaxial  shaker 
was  presented  previously^. 

For  the  problem  of  the  triaxial  shaker  system,  a  model 
was  constructed  that  represented  the  table,  the  eight  ac¬ 
tuators,  and  the  eight  bearings.  Nominal  values  for  the 
coefficients  in  this  model  were  derived  from  a  finite  el¬ 
ement  model  and  from  manufacturer’s  data.  Unfortu¬ 
nately,  the  model  was  only  moderately  predictive  of  the 
actual  response  of  the  system  over  a  300  Hz  band  for 
which  test  data  was  available.  Three  approaches  were 
evaluated  for  the  purpose  of  updating  the  model  coeffi¬ 
cients  in  order  to  obtain  better  agreement  with  the  test 
data.  Two  are  based  on  input  (or  equation)  error  formu¬ 
lations,  while  the  third  is  based  on  an  output  error  for¬ 
mulation.  These  approaches  are  discussed  in  section  3. 
They  have  been  applied  to  two  different  test  cases;  the 
results  are  discussed.  Work  continues  on  applying  these 
approaches  to  the  actual  system. 

A  second  modeling  approach  was  evaluated  in  an  at¬ 
tempt  to  circumvent  problems  encountered  using  a 
physically-based  model.  This  second  approach  involves 
fitting  the  system  response  with  a  parametric  model  of 
the  system,  in  this  case,  an  autoregressive  moving  aver¬ 
age  (ARMA)  model. 

Many  methods  for  parametric  modeling  have  been  pro¬ 
posed.  Several  involve  time-domain  measurements.  For 
instance,  Hu,  Chen,  and  Wu  considered  the  case  where 
a  tiine-domain  input/output  pair,  or  just  an  output,  is 
available  for  every  degree  of  freedom  of  the  system^ 
They  used  a  modified  Yule- Walker  algorithm  to  obtain 
autoregressive  coefficients,  which  were  then  used  to  pro¬ 
duce  modal  parameters.  Jabbari  and  Gibson  used  a 
lattice  filter  to  estimate  the  order  and  coefficients  of 
an  ARMA  model  given  time-domain  input  and  output 
measurements^^ 

If  the  available  measurements  are  in  the  frequency  do¬ 
main  (this  is  the  case  for  the  much  of  the  available  triax¬ 
ial  shaker  system  data),  other  approaches  are  available. 
Chen,  Juang,  and  Lee  converted  frequency  response  data 
to  Markov  parameters  by  way  of  a  finite- ordered  matrix 


fraction,  and  then  used  the  Markov  parameters  to  con¬ 
struct  a  linear  state  space  model  by  way  of  a  eigensys- 
tem  realization  algorithm  (which  uses  a  singular  value 
decomposition)^^.  Jean-Christophe  obtained  a  transfer 
function  from  modal  data^^.  Again,  the  literature  con¬ 
tains  many  other  possible  approaches. 

For  this  problem,  the  approach  of  Friedlander  and 
20 

For  at  proved  useful.  They  proposed  using  a  modified 
Yule- Walker  algorithm  to  generate  ARMA  coefficients 
(for  an  assumed  order)  from  frequency  response  data. 
The  entire  algorithm  is  quite  compact  computationally; 
it  involves  the  least-squares  solution  of  three  linear  sys¬ 
tems  and  a  fast  Fourier  transform.  The  technique  is 
briefly  described  in  section  4.  The  results  of  deriving 
models  using  this  technique  are  also  presented. 

Conclusions  regarding  the  simulation  of  the  triaxial 
shaker  system  and  the  techniques  used  to  derive  the 
models  are  presented  in  section  5. 


2.  Simulation  Design 

An  engineering  simulation  of  the  triaxial  shaker  system  is 
being  constructed.  The  system’s  actuators  are  Ling  Elec¬ 
tronics  model  2106  electrodynamic  shakers  connected  to 
double-spherical  hydrostatic  bearings,  which  are  in  turn 
connected  to  the  40-inch  square  magnesium  table.  Each 
shaker  consists  of  a  coil,  mounted  on  a  fixed  base,  and 
an  armature. 

Two  different  models  have  been  implemented.  The  first 
is  physically-based  model  consisting  of  a  table  model  and 
models  for  the  eight  actuators,  including  the  bearings. 
The  table  is  modeled  by  a  linear  system  of  springs  and 
masses  derived  from  a  finite  element  model.  The  origi¬ 
nal  form  of  this  system  contained  7,392  structural  nodes 
for  a  total  of  44,352  degrees  of  freedom^^  Since  this 
results  in  mass  and  stiffness  matrices  with  nearly  two  bil¬ 
lion  elements  each,  this  system  is  not  practical  for  com¬ 
puter  simulation.  Instead,  the  table  system  was  reduced 
to  66  degrees  of  freedom  through  a  generalized  dynamic 
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reduction  (an  extension  to  Guyan  reduction)  The 

66  degrees  of  freedom  represent  six  degrees  of  freedom  at 
each  of  the  eight  actuator  attachment  points  (known  as 
control  points)  and  eighteen  additional  degrees  of  free¬ 
dom  used  to  fit  the  original  modal  data.  The  reduced 
model  matches  the  full  model  to  within  one  percent  in 
frequency  for  each  of  the  first  four  bending  modes  and  to 
within  five  percent  for  each  of  the  next  six  modes.  The 
models  for  the  actuators  consist  of  a  system  of  springs, 
masses,  and  dampers.  Values  for  the  various  coefficients 
were  estimated  using  manufacturer’s  data  and  the  tech¬ 
niques  described  in  section  3. 

The  table  bending  model,  even  after  reduction,  is  large, 
has  many  high  frequency  modes,  and  is  computationally 
expensive  to  integrate.  Rather  than  simulate  all  of  the 
modes,  only  a  subset  is  actually  simulated.  First,  the  ta¬ 
ble  coordinates  are  transformed  into  modal  coordinates. 
Only  the  modal  states  with  natural  frequencies  of  less 
than  some  cutoff  frequency  are  actually  integrated;  the 
others  are  set  to  zero.  The  cutoff  frequency  is  chosen  as 
five  times  the  highest  frequency  of  interest. 

A  fixed-step-size  Runge-Kutta  scheme  was  used  as  the 
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baseline  numerical  integrator  Even  after  the  highest 
frequency  bending  modes  were  neglected,  this  method 
required  a  small  step  size.  Experiments  were  conducted 
with  other  schemes,  including  a  variable-step  Runge- 
Kutta  scheme,  the  Adams-Moulton  method,  and  Gear’s 
backward  differentiation  formula  (BDF)  method^^ 

For  low  frequency  inputs,  the  high  frequency  modes  de¬ 
cay  after  only  a  fraction  of  a  period  and  the  system 
becomes  numerically  stiff  .  The  implicit  methods 
(Adams-Moulton  and  Gear’s  BDF)  tended  to  be  more  ef¬ 
ficient  after  the  initial  transients  decayed,  but  displayed 
a  sensitivity  to  sudden  changes  in  the  external  forces. 

The  second  model  has  the  form  of  an  autoregressive  mov¬ 
ing  average  function.  The  coefficients  were  derived  from 
available  test  data,  as  described  in  section  4.  Since  the 
ARMA  model  produces  the  system  response  directly,  no 
numerical  integration  is  necessary,  and  the  simulation 
executes  rapidly. 

A  simple  control  algorithm  was  implemented  in  the  simu¬ 
lation.  This  allowed  for  pfeliniinary  studies,  in  particular 
comparisons  to  actual  testing,  to  be  conducted.  Even¬ 
tually,  the  actual  control  algorithms  used  in  the  triaxial 
shaker  system  will  be  integrated  into  the  simulation. 

3.  Identification  of  a  Physically-Based  Model 

An  initial  physically-based  model  of  the  triaxial  shaker 
system  was  developed  from  a  finite  element  analysis  and 
manufacturer  data.  After  linearization,  the  model  as¬ 
sumed  a  system  of  the  form 

X  =  A(0)x -j- B(0)u  (1) 

y  =  C(0)x  +  B(0)u  -I  V  (2) 


where  x  is  the  state,  u  is  the  control,  y  is  the  measure¬ 
ment,  and  V  is  measurement  noise.  In  particular  for  this 
work,  the  state  vector  x  consists  of  positions  and  ve¬ 
locities,  and  the  measurements  are  accelerations.  The 
matrices  A,  B,  C,  and  D  depend  on  the  vector  of  pa¬ 
rameters  6  that  are  to  be  updated.  For  this  work,  the 
matrix  D  is  assumed  to  be  zero. 

Three  algorithms  for  estimating  9  are  considered.  The 
first  uses  an  output  error  formulation.  Let  Gik{uJn)  be 
the  measured  frequency  response  of  the  system  from  the 
kth  input  control  to  the  ith  output  acceleration,  at  AT  — 
1  frequencies  (jJn-  The  model  frequency  response  from 
input  k  to  output  i  is 

Gikii^n)  =  Ci{0)  {ju>ni  -  A{e)y^  Bk{e)  (3) 

The  cost  function  for  the  first  algorithm  is 

N—  Ip  q 

n=0  i=l  A:=l 

where  q  is  the  number  of  inputs  and  p  is  the  num¬ 
ber  of  outputs.  The  cost  function  is  minimized  by  a 
Nelder-Mead  simplex  method,  which  requires  no  gradi¬ 
ent  information^^. 

The  second  and  third  algorithms  use  input  error  formu¬ 
lations.  Let 

Sk{u^„,e)  =  {juni  -  A{e))  -  Bk{e)  (s) 

for  1  <  /c  <  g',  where  dX/dUk  is  estimated  from  the 
measurements^.  Then  the  cost  function  for  the  second 
algorithm  is 

N-l  q 

(6) 

n=0  fc=l 

The  cost  function  J2{0)  is  minimized  with  respect  to  9 
through  a  Nelder-Mead  simplex  search. 

The  third  algorithm  uses  the  same  cost  function,  J2(^), 
as  the  second  algorithm,  but  minimizes  it  through  a 
least-squares  approach^. 

The  three  algorithms  were  applied  to  two  different  test 
cases  to  evaluate  their  relative  performance;  those  results 
follow.  Work  continues  in  applying  the  most  promising 
of  these  techniques  to  actual  shaker  test  data. 

3.1.  The  First  Test  Case 

The  physical  system  for  the  first  test  case  is  depicted 
in  figure  2.  Each  small  block  represents  a  spring  and 
damper  in  parallel.  The  stiffness  and  damping  matrices, 
Ki  and  Ci,  are  composed  of  linear  combinations  of  the 
spring  and  damper  coefficients  for  the  blocks  indicated 
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in  the  figure.  The  mass  matrix,  Mi,  is  a  diagonal  matrix. 
The  system  matrices  are  then 


A 

B 

C 


0  I 

0000  Pi  000 

0000  0  P2OO 

A7i 


(7) 

(8) 
(9) 


where  pi  and  p2  are  control  gains  and  where  Aji  and  Agi 
represent  the  seventh  and  eighth  rows  of  A. 


parameter 

percent  error  versus  truth  || 

initial 

OE/S 

lE/S 

lE/LS 

if  1  (1,1) 

15.0 

-8.7 

-6.4 

0.0 

^1(1,  2) 

15.0 

12.8 

23.1 

0.0 

Ki(l,3) 

15.0 

23.9 

42.3 

0.0 

15.0 

10.7 

53.0 

0.0 

if  1(2,  2) 

15.0 

18.7 

9.2 

0.0 

Ki(2,3) 

15.0 

-28.3 

-2.3 

0.0 

Kj{2,  4) 

15.0 

-36.7 

2.2 

0.0 

i^l(3,3) 

15.0 

-17.2 

3.1 

0.0 

Ki(3,4) 

15.0 

11.3 

-3.2 

0.0 

i^l(4,4) 

15.0 

3.7 

2.6 

0.0 

Ci(l,l) 

-15.0 

34.4 

17.2 

0.0 

Ci(l,2) 

-15.0 

13.8 

11.8 

0.0 

Ci(l,3) 

-15.0 

-29.3 

-10.7 

0.0 

Cl  (1,4) 

-15.0 

-21.2 

15.2 

0.0 

Cl  (2,  2) 

-15.0 

17.0 

-10.5 

0.0 

Cl  (2,  3) 

-15.0 

0.6 

-2.6 

0.0 

^1(2,4) 

1  -15.0 

-20.0 

11.2 

0.0 

Cl  (3, 3) 

-15.0 

23.0 

12.9 

0.0 

Cl  (3,  4) 

-15.0 

-20.7 

-133.1 

0.0 

Cl  (4,  4) 

-15.0 

-1.5 

7.1 

0.0 

91 

66.7 

-6.2 

-26.7 

0.0 

92 

100.0 

46.0 

-0.1 

0.0 

Table  1;  Estimated  Parameters 


squares  solution  to  the  input  error  formulation  produces 
gorul  estimates  of  the  unknown  parameters.  The  esti- 
mafos  of  the  parameters  using  the  simplex  solutions  to 
hf.ih  the  input  error  and  output  error  formulations  are 
not  as  good,  though  the  frequency  response  produced 
by  the  erroneous  models  is  still  close  to  the  measured 
re-';ponse.  This  suggests  that  the  search  algorithms  con¬ 
verged  on  local  minima  away  from  the  actual  parameter 
values. 

Tlie.  least  squares  solution  to  the  input  error  formulation 
(IK/'LS)  converged  in  less  than  ten  iterations.  The  con¬ 
vergence  for  the  simplex  minimization  of  both  the  input 
(If’./S)  and  output  (OE/S)  formulations  was  somewhat 
slower,  each  requiring  several  thousand  evaluations  of  the 
cont  function  (figure  4).  The  output  error  formulation 
converged  more  slowly  than  the  input  error  formulation. 


Figure  3:  Test  Case  1  Model  Response 


Results  for  the  three  algorithms  for  the  first  test  case  are 
presented  in  table  1  and  figures  3  and  4.  The  parameters 
chosen  for  this  test  case  were  the  elements  of  the  symmet¬ 
ric  matrices  Ki  and  Ci  and  the  control  gains.  The  least 


Figure  4:  Test  Case  1  Convergence 
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3.2.  The  Second  Test  Case 


/rpresont  the  sixth,  seventh,  and  eighth  rows  of  A. 


The  second  test  case  is  a  portion  of  a  model  for  a  triaxial 
shaker  system.  The  actual  system  consists  of  eight  actu¬ 
ators  attached  to  a  40  inch  by  40  inch  square  table.  Four 
of  the  actuators  are  attached  to  the  bottom  of  the  table 
and  are  aligned  vertically,  while  the  other  four  actuators 
are  attached  to  the  sides  and  are  aligned  horizontally. 
Because  the  immediate  objective  is  to  evaluate  the  rela¬ 
tive  performance  of  the  three  identification  algorithms,  a 
simplified  model  of  only  one  axis  of  the  horizontal  motion 
is  considered.  Note,  however,  that  it  is  straightforward. 

n 

to  extend  the  model  to  three  dimensions  . 


The  physical  model  for  the  second  test  case  is  depicted 
in  figure  5.  Each  small  block  represents  a  spring  and 
damper  in  parallel.  Only  forces  and  motion  along  the 
axis  indicated  by  the  arrows  are  considered.  The  springs 
and  dampers  are  assumed  linear  in  that  direction.  The 
masses  represent  the  mass  of  the  shaker  table,  while  the 
springs  and  dampers  represent  table  and  actuator  com¬ 
pliance  and  damping.  The  gain  is  the  actuator  gain,  and 
the  measurements  are  from  accelerometers  positioned  as 
indicated.  The  stiffness  and  damping  matrices,  K2  and 
C2,  are  composed  of  linear  combinations  of  the  individ¬ 
ual  spring  and  damper  coefficients,  ki  and  Ci,  for  each 
block  indicated  in  the  figure.  The  mass  matrix,  M2,  is  a 
diagonal  matrix. 


Figure  5:  Test  Case  2  Schematic 


The  system  matrices  are  then 


A 

B 

C 


0  / 
~M-^K2  -M^^C2 

OOOOi^OOO 


Aqi 

Aji 

A%i 


(10) 

(11) 

(12) 


where  p  is  a  control  gain  and  where  Aei, 


Avi,  and  Asi 


fiaranieter 

percent  error  versus  truth  {| 

initial 

OE/S 

lE/S 

lE/LS 

^'1 

15.00 

-0.01 

0.01 

0.00 

A'2 

15.00 

-0.07 

0.02 

0.00 

15.00 

0.67 

-0.17 

0.00 

/C4 

15.00 

0.02 

-0.00 

0.00 

Cl 

-15.00 

-0.19 

-0.03 

0.00 

C2 

-15.00 

-0.52 

0.33 

0.00 

C3 

-15.00 

0.75 

0.18 

0.00 

C4 

-15.00 

-0.22 

-0.06 

0.00 

9 

233.33 

-0.08 

0.00 

0.00 

Table  2:  Estimated  Parameters 


Results  for  the  three  algorithms  for  the  second  test  case 
are  presented  in  table  2  and  figures  6  and  7.  The  param- 
^:ters  chosen  for  this  test  case  were  the  spring  coefficients, 
the  damping  coefficients,  and  the  control  gain.  All  three 
algorithms  produce  good  estimates  of  the  unknown  pa- 
ra.ineters.  The  resulting  frequency  responses  are  likewise 
good. 


The  lecist  squares  solution  to  the  input  error  formulation 
converged  in  less  than  ten  iterations.  The  convergence 
for  the  simplex  minimization  of  the  input  error  formu¬ 
lation  took  hundreds  of  evaluations  of  the  cost  function, 
while  the  output  error  formulation  required  several  thou¬ 
sand  (figure  7). 

A  comparison  of  the  computational  effort  required  for 
eacli  algorithm,  for  both  test  cases,  is  given  in  figure  8. 
The  computational  effort  for  the  output  error  algorithm 
is  driven  by  the  need  to  repeatedly  invert  the  matrix 
in  (3),  while  the  effort  for  the  two  input  error  algorithms 
is  driven  by  the  need  to  pseudoinvert  another  matrix. 
Although  the  simplex  minimization  of  the  input  error 
formulation  requires  fewer  cost  function  evaluations  than 
the  output  error  approach,  the  expense  of  the  pseudoin¬ 
version  results  in  a  greater  overall  computational  burden. 
Because  the  least-squares  minimization  of  the  input  er¬ 
ror  formulation  converges  in  so  few  steps,  its  total  cost 
is  lowest,  in  spite  of  the  pseudoinversion.  Note  also  that 
the.  dependence  on  a  pseudoinversion  makes  the  input 
error  formulations  expensive  for  system  matrices  with  a 
large  dimension. 

Test  case  2  required  less  computational  effort  than  test 
case  1.  This  is  due  in  part  to  the  parameterization  of 
the  problems.  For  test  case  1,  no  assumption  was  made 
regarding  the  structure  of  the  stiffness  and  damping  ma¬ 
trices,  other  than  that  of  dimension.  The  elements  them- 
seh'es  of  these  matrices  were  taken  as  the  unknown  pa¬ 
rameters.  For  test  case  2,  a  certain  structure,  including 
some  symmetry,  was  assumed  for  the  stiffness  and  damp¬ 
ing  matrices,  and  the  parameters  included  unknown  val¬ 
ues  within  that  structure.  This  not  only  decreased  the 
number  of  unknown  parameters  (relative  to  test  case  1), 
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but  it  also  constrained  the  search  for  a  minimizing  pa¬ 
rameter  set  in  such  a  way  as  to  encourage  convergence. 

The  simplex  method  of  minimizing  a  cost  function  is 
probably  not  the  most  efficient  algorithm,  but  it  has 
the  advantage  of  simplicity  in  implementation.  For 
this  problem,  it  also  appeared  to  have  the  advantage 
of  being  better  able  to  avoid  local  minima,  as  com¬ 
pared  to  other  techniques  that  were  considered,  includ¬ 
ing  a  Broyden-Fletcher-Goldfarb-Shanno  variable  metric 
(quasi-Newton)  method 
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frequency  (Hz) 


Figure  6:  Test  Case  2  Model  Response 


evaluation  number 


Figure  7:  Test  Case  2  Convergence 

4.  Identification  of  a  Parametric  Model 

Rather  than  using  a  model  based  on  a  finite  element 
model  of  the  table  and  models  for  the  actuators,  an  au¬ 
toregressive  moving  average  (ARMA)  model  can  be  used 
for  each  of  the  input/output  pairs.  This  approach  results 
in  a  model  that  compares  well  with  the  test  data.  The 
disadvantages  of  this  approach  are  that  the  parameters 
themselves  have  no  physical  significance  and  that  there 
is  no  convenient  way  to  incorporate  a  priori  knowledge 
(specifically,  the  finite  element  model)  into  the  model. 
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Method 


Figure  8:  Computational  Effort 


In  this  approach,  the  identification  problem  is  to  find  a 
model  transfer  function  G(^),  which  has  the  form 


G{z) 

^  6^ 
a{z) 

P 

(13) 

a{z) 

—  aiZ~^ 

i=l 

P 

(14) 

h{z) 

i=l 

(15) 

such  that  the  error 

e 

—  Gk  —  Gk 

(16) 

is  minimized,  where 

Gk  = 

(17) 

Gk  = 

(18) 

for  fc  =  0, 1, . . . ,  A/"  —  1.  The  true  transfer  function  G{z) 
is  assumed  to  be  of  the  form 

G{z) 

II 

(19) 

a{z) 

Q 

(20) 

b[z) 

T 

II 

(21) 

on  ^ 

The  procedure  used  is  that  of  Friedlander  First,  an 
unwindowed  correlation  sequence  is  generated  by  taking 
the  discrete  inverse  Fourier  transform  of  the  square  of 
the  desired  amplitude  response.  Next,  a  windowed  cor¬ 
relation  sequence  is  formed  by  multiplying  the  unwin¬ 
dowed  sequence  by  a  Hamming  window,  and  then  the 
autoregressive  coefficients  a  are  computed.  In  order  to 
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compute  the  moving  average  coefficients  6,  an  additive 
decomposition  is  performed,  followed  by  the  construction 
of  the  complex  cepstrum  of  the  estimated  transfer  func- 

on 

tion  G{z)  ,  from  which  the  moving  average  coefficients 
are  derived. 

Figures  9  through  12  show  a  comparison  between  test 
data  and  the  ARM  A  model.  Figures  9  and  10  compare 
the  frequency  response  from  test  data  and  the  model 
frequency  response  for  three  cases,  for  p  equal  to  5,  10, 
and  20,  for  testing  done  with  an  empty  table.  The  model 
shows  good  agreement  for  p  =  5  and  excellent  agreement 
for  higher  p.  Results  for  the  other  fourteen  input/output 
pairs  were  similarly  good.  (Due  to  symmetry,  the  re¬ 
sponses  of  the  eight  control  points  to  the  eight  actuators 
can  be  characterized  by  sixteen  input/output  pairs:  eight 
responses  to  a  horizontal  actuator  and  eight  responses  to 
a  vertical  actuator.) 

Figures  11  and  12  show  the  results  for  a  table  with  a 
payload  attached.  The  payload  is  a  test  fixture  for  an 
NS 20  guidance  system.  Once  again,  the  model  shows 
good  agreement  with  the  test  data. 


Figure  9:  ARM  A  Estimation,  Empty  Table  (Vertical  Ex¬ 
citation) 

5.  Conclusions 

An  engineering  simulation  of  the  triaxial  shaker  system 
allows  for  the  evaluation  of  system  capabilities  and  limi¬ 
tations,  the  support  of  test  methodology  and  test  proce¬ 
dure  development,  and  the  testing  of  controller  software 
modifications.  The  implementation  of  the  system  models 
is  designed  to  produce  rapid  execution. 

Two  different  modeling  approaches  were  considered.  For 
the  first  approach,  which  used  a  physically-based  model, 
two  different  formulations  were  used  for  updating  the 
model  with  respect  to  test  data.  The  output  error  for¬ 
mulation  is  hindered  by  the  presence  of  significant  non- 


Figure  10:  ARMA  Estimation,  Empty  Table  (Horizontal 
Excitation) 


linearities  in  its  cost  function,  and  it  displayed  only  mod¬ 
erate  performance  for  the  test  cases.  The  input  error  for¬ 
mulations  performed  very  well  for  the  two  test  cases.  A 
least-squares  approach  to  minimizing  the  cost  function 
produced  by  the  input  error  formulation  was  particularly 
effective.  Work  continues  on  applying  the  most  promis¬ 
ing  of  these  techniques  to  the  triaxial  shaker  system. 

The  second  modeling  approach  used  a  parametric  model, 
in  particular,  an  autoregressive  moving  average  model. 
This  technique  produced  a  model  that  agreed  quite  well 
vith  the  test  data.  The  major  disadvantage  to  the 
parametric  approach  is  that,  unlike  the  physically-based 
model,  there  is  no  direct  physical  interpretation  of  the 
resulting  coefficients.  This  precludes  using  the  paramet¬ 
ric  model  for  studies  of  the  system  behavior  under  per¬ 
turbations  of  physical  characteristics.  Nevertheless,  this 
approach  is  useful  for  evaluating  the  capabilities  and  lim¬ 
itations  of  the  triaxial  shaker  system. 
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